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The double pulsar PSR J0737— 3039A/B consists of two neutron stars in a 
highly relativistic orbit that displays a roughly 30-second eclipse when pulsar 
A passes behind pulsar B. Describing this eclipse of pulsar A as due to absorp- 
tion occurring in the magnetosphere of pulsar B, we successfully use a simple 
geometric model to characterize the observed changing eclipse morphology 
and to measure the relativistic precession of pulsar B's spin axis around the 
total orbital angular momentum. This provides a test of general relativity 
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and alternative theories of gravity in the strong-field regime. Our measured 
relativistic spin precession rate of 4.77to!65 °J T _1 (68% confidence level) is con- 
sistent with that predicted by general relativity within an uncertainty of 13%. 

Spin is a fundamental property of most astrophysical bodies, making the study of its grav- 
itational interaction an important challenge (1). Spin interaction manifests itself in different 
forms. For instance, we expect the spin of a compact rotating body in a binary system with 
another compact companion to couple gravitationally with the orbital angular momentum (rel- 
ativistic spin-orbit coupling) and also with the spin of this companion (relativistic spin-spin 
coupling) (|2]|21)- Observing such phenomena provides important tests for theories of gravity, 
because every successful theory must be able to describe the couplings and to predict their 
observational consequences. In a binary system consisting of compact objects such neutron 
stars, one can generally consider the spin-orbit contribution acting on each body to dominate 
greatly the spin-spin contribution. This interaction results in a precession of the bodies' spin 
axis around the orbital angular momentum of the system, behavior we refer to as relativistic 
spin precession. 

While relativistic spin precession is well studied theoretically in general relativity (GR), 
the same is not true of alternative theories of gravity and hence, quantitative predictions of 
deviations from GR spin precession do not yet exist (2). For instance, it is expected that in 
alternative theories relativistic spin precession may also depends on strong self-gravitational 
effects, i.e. the actual precession may depend on the structure of a gravitating body (2). In the 
weak gravitational fields encountered in the solar system, these strong-field effects generally 
cannot be detected (|5]|6][7|). Measurements in the strong-field regime near massive and compact 
bodies such as neutron stars and black holes are required. Relativistic spin precession has been 
observed in some binary pulsars [e.g. (|^ll9l [70|) l. but it has usually only provided a qualitative 
confirmation of the effect. Recently, the binary pulsar PSR B 1534+ 12 has allowed the first 
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quantitative measurement of this effect in a strong field, and although the spin precession rate 
was measured to low precision, it was consistent with the predictions of GR (11). 

Here we report a precision measurement of relativistic spin precession using eclipses ob- 
served in the double pulsar C/2]|75]). This measurement, combined with observational access to 
both pulsar orbits in this system, allows us to constrain quantitatively relativistic spin precession 
in the strong-field regime within a general class of gravitational theories that includes GR. 

PSR J0737— 3039A/B consists of two neutron stars, both visible as radio pulsars, in a rela- 
tivistic 2.45-hour orbit (|72II73I) . High-precision timing of the pulsars, having spin periods of 23 
ms and 2.8 s (hereafter called pulsars A and B, respectively), has already proven to be the most 
stringent test-bed for GR in the strong-field regime ([7?]). and enables four independent timing 
tests of gravity, more than any other binary system. 

The orbital inclination of the double pulsar system is such that we observe the system al- 
most perfectly edge-on. This coincidence causes pulsar A to be eclipsed by pulsar B at pulsar 
A's superior conjunction (13). The modestly frequency-dependent eclipse duration, about 30 
s, corresponds to a region extending ~1.5x 10 7 m (|75l) . The light curve of pulsar A during its 
eclipse shows flux modulations that are spaced by half or integer numbers of pulsar B's rota- 
tional period (|7oT) . This indicates that the material responsible for the eclipse corotates with 
pulsar B. The relative orbital motions of the two pulsars and the rotation of pulsar B thus allow 
a probe of different regions of pulsar B's magnetosphere in a plane containing the line of sight 
and the orbital motion. 

Synchrotron resonance with relativistic electrons is the most likely mechanism for efficient 
absorption of radio emission over a wide range of frequencies. In the model proposed by Lyu- 
tikov and Thompson (1 ), this absorbing plasma corotates with pulsar B and is confined within 
the closed field lines of a magnetic dipole truncated by the relativistic wind of pulsar A. The 
dipole magnetic moment vector makes an angle a with respect to the spin axis of pulsar B, 



whose orientation in space can be described by two angles: the colatitude of the spin axis with 
respect to the total angular momentum of the system, 9, and the longitude of the spin axis, <p 
(see Fig. 1 for an illustration of the system geometry). Additional parameters characterizing the 
plasma opacity, /x, the truncation radius of the magnetosphere, -R mag , and the relative position 
of pulsar A with respect to the projected magnetosphere of pulsar B, z , are also included in the 
model ©. 

We monitored the double pulsar from December 2003 to November 2007 using the Green 
Bank Telescope in West Virginia; most of the data were acquired as part of the timing observa- 
tions reported in (14). The data used for our analysis were taken at 820 MHz with the SPIGOT 
instrument (18), which provides 1024 frequency channels across a 50 MHz bandwidth. Data 
for a total of 63 eclipses of pulsar A were collected over the 4-year period, with many ob- 
tained during semi-annual concentrated observing campaigns. We dedispersed each eclipse 
data set by adding time shifts to frequency channels in order to compensate for the frequency- 
dependent travel time of radio waves in the ionized interstellar medium and we then folded 
them at the predicted spin period of pulsar A using the pulsar analysis packages PRESTO (|79|) 
and SIGPROC (20) (see (14) for details about the radio timing). Next, we extracted the rela- 
tive pulsed flux density of pulsar A by fitting each folded interval for the amplitude of a high 
signal-to-noise ratio pulse profile template made from the integrated pulse observed during the 
several-hour observation that includes each eclipse. Finally, we normalized the flux densities so 
the average level outside the eclipse region corresponded to unity. We chose the time resolution 
of our eclipse light curves to equal, on average, four individual pulses of pulsar A (~91 ms). 

In addition to the flux density, we determined the orbital phase and the spin phase of pulsar 
B corresponding to each data point of our time series. Orbital phases were derived from the 
ephemeris published in (14} . Spin phases were empirically measured from data folded at the 
predicted period of pulsar B in a way similar to that described above for pulsar A. Over the four- 



year monitoring campaign, we found significant changes in pulsar B's pulse profile, likely due 
to the precession of its spin axis, which were also reported in (|2i|) . Around 2003, the average 
pulse profile was unimodal, resembling a Gaussian function. It evolved such that by 2007, it 
displayed two narrow peaks. Using the pulse peak maximum as a fiducial reference point is 
certainly not appropriate. We find, however, that the unimodal profile gradually became wider 
and then started to form a gap near the center of its peak. Since then, the outer edges of the 
pulse profile have not significantly changed but the gap evolved such that two peaks are now 
visible. This lets us presume that the underlying average profile is reminiscent of a Gaussian- 
like profile to which some "absorption" feature has been superimposed near the center, leaving 
a narrow peak on each side. We therefore defined the fiducial reference point to lie at the center 
of the unimodal "envelope" that we reconstructed from the first ten Fourier bins of the pulse 
profile, which contains 512 bins in total (see Fig. S2 of the Supporting Online Material for an 
illustration of the pulse profile evolution). 

We implemented the eclipse modeling of our data in two steps: the fitting of individual 
eclipse profiles and the search for evolution of the geometry of pulsar B. We first searched the 
full phase space to identify best-fit values of six parameters (see Supporting Online Material for 
more details). Then, we reduced the number of free parameters to the subset (9, <fi, a) describing 
the orientation of pulsar B's spin and magnetic axes by fixing the other parameters to their best- 
fit values: [i = 2, -R mag = 1.29° (projected value in terms of orbital phase) and z /R mag = 
—0.543 (see Fig. 1). Finally, we performed a high-resolution mapping of the likelihood of 
this subspace in order to investigate subtle changes in the geometry. Lyutikov and Thompson 
© predicted that such changes, due to relativistic spin precession, could affect the eclipse 
light curve. In principle, relativistic spin precession of pulsar B's spin axis around the total 
angular momentum should induce a secular change of the longitude of the spin axis, <p, while 
the magnetic inclination, a, and the colatitude of the spin axis, 9, are expected to remain fixed 



over time. Indeed, from model fitting, we find no significant time evolution of a and 6, whereas 
(ft does change. Because of correlation between the parameters, we jointly evaluated the best-fit 
geometry of pulsar B using a time-dependent model in which a = a and 9 = 9 are constants, 
and (p varies linearly with time, i.e. <\> = (fi — VL B t, where f2 B is the rate of change of pulsar 
B's spin axis longitude and the epoch of 4> = (fio is May 2, 2006 (MJD 53857). Figure 2 
shows the time evolution of the parameters and the fit derived from this joint time-dependent 
model (Table [I])- The precession rate f2 B of 4.77^0^5 °yr _1 (|22|) agrees with the precession 
rate predicted by GR (EH), 5.0734 ± 0.0007 °yr _1 d23), within an uncertainty of 13% (68% 
confidence level). 

This relatively simple model (pO) is able to reproduce the complex phenomenology of the 
eclipses (see Fig. 3 and Movie SI in Supporting Online Material) except at the eclipse bound- 
aries where slight magneto spheric distortions or variations in plasma density are likely to occur. 
Fits including the egress generally are poor in the central region where we observe narrow 
modulation features, which are critical for determining pulsar B's geometry. For this reason, 
we excluded the egress from the fits, using orbital phases between —1.0 and 0.75° (see Fig. 3). 
We accounted for systematics introduced by the choice of the region to fit in the priors of our 
Bayesian model (see Supporting Online Material). This improved the fit of the model through- 
out the center region of the eclipse while still producing qualitatively good predictions near the 
eclipse egress. The overall success of the model implies that the geometry of pulsar B's magne- 
tosphere is accurately described as predominantly dipolar; a pure quadrupole, for instance, does 
not reproduce the observed light curves. Although the model does not exclude the possibility 
that higher-order multipole components may exist close to the surface of pulsar B, our model- 
ing supports the conclusions (1) that these eclipses yield direct empirical evidence supporting 
the long-standing assumption that pulsars have mainly dipolar magnetic fields far from their 
surface. 



The direct outcome from modeling the eclipse profile evolution is a measurement of the 
effect of relativistic spin precession (see Movie S2 in Supporting Online Material for an il- 
lustration of the time evolution of the eclipse). We can use the inferred precession rate to 
test GR (see Fig. 4) and to further constrain alternative theories of gravity and the strong- 
field aspects of relativistic spin precession. We use the generic class of relativistic theories 
that are fully conservative (Lorentz-invariant) and based on a Lagrangian, as introduced by 
Damour & Taylor ©. In this way we can study the constraints of our observations on the- 
ories of gravity by describing the spin-orbit interaction within a specific theory by coupling 
functions appearing in the corresponding part of the Lagrangian. In this framework, we can 
write the precession rate of pulsar B in a general form, f2e = ctbL/ ofj(l — e 2 ) 3 / 2 where L 
is the orbital angular momentum of the system, a B is the semimajor axis of the relative orbit 
between the pulsars, e the eccentricity of the orbit and a B is a generic strong-field spin-orbit 
coupling constant. Since L and <ir are not directly measurable, it is more convenient to write 
the above expression using observable Keplerian and post-Keplerian parameters. While alter- 
native forms generally involve a mixture of gravitational theory-dependent terms, the particular 
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choice Vt-Q = XA s 2 B x -^-^ x ^-^p- is the only one that does not incorporate further theoretical 
terms other than the spin-orbit coupling constant, a B , the speed of light, c, and a generalized 
gravitational constant for the interaction between the two pulsars, Q. In this expression, the 
Keplerian parameters e and n = 2n/Pb, the angular orbital frequency, are easily measurable 
for any binary system. On the other hand, the post-Keplerian Shapiro delay shape parameter s, 
equivalent to the sine of the orbital inclination angle (2), requires relatively edge-on orbits to 
be observed. Measurement of the projected semi-major axes of the two orbits (|25l) . xa and x B , 
found in the above equation, necessitates that each body must be timeable. Therefore, the dou- 
ble pulsar is the only relativistic binary system that allows a direct constraint on the spin-orbit 
coupling in general theories of gravity. Using the inferred precession rate of fie = 4.77^a65 ° 



yr _1 , we derive ( £ -§ E -) = 3.38ig;4g. Every successful theory of gravity in the given generic 
framework must predict this value — these observations provide a strong-field test of gravity 
that complements and goes beyond the weak-field tests of relativistic spin precession (1261) . In 
GR, we expect to measure (^p) = 2 + l^t = 3.60677 ± 0.00035, where we have used 
the masses determined from the precisely observed orbital precession and the Shapiro delay 
shape parameter under the assumption that GR is correct (|7?1) . Comparing the observed value 
with GR's predictions, we find (^p) / (^p) = 0.94 ± 0.13. Hence, GR passes this test 
of relativistic spin precession in a strong-field regime, confirming, within uncertainties, GR's 
effacement property of gravity even for spinning bodies, i.e. the notion that strong internal grav- 
itational fields do not prevent a compact rotating body from behaving just like a spinning test 
particle in an external weak field (1271) . 

The spin precession rate, as well as the timing parameters entering in the calculation of 
( ^iP" ) ' are a ^ ^dependent of the assumed theory of gravity. If the main contribution limiting 
the precision of this new strong-field test comes from the inferred spin precession rate, we ex- 
pect that the statistical uncertainty should decrease significantly with time, roughly as the square 
of the monitoring baseline for similar quantity and quality of eclipse data. The contribution of 
systematics to the error budget should also decrease, but its functional time dependence is dif- 
ficult to estimate. Although the orbital and spin phases of pulsar B are input variables to the 
eclipse model, our ability to determine the orientation of pulsar B in space does not require the 
degree of high-precision timing needed for measurement of post-Keplerian parameters; evalu- 
ating spin phases to the percent-level, for instance, is sufficient. Therefore, the intrinsic cor- 
rectness of the model and its ability to reproduce future changes in the eclipse profile due to 
evolution of the geometry are the most likely limitations to improving the quality of this test of 
gravity, at least until the measured precession rate reaches a precision comparable with the tim- 
ing parameters involved in the calculation of \jg-\ Better eclipse modeling could be achieved 



from more sensitive observations and thus new generation radio telescopes such as the proposed 
Square Kilometer Array could help make important progress. Pulsar A does not show evidence 
of precession (T2ffll29l) likely because its spin axis is aligned with the orbital angular momentum; 
it should therefore always remain visible, thus allowing long-term monitoring of its eclipses. 
Pulsar B, however, could disappear if spin precession causes its radio beam to miss our line of 
sight (|271) . In this event, we would need to find a way to circumvent the lack of observable spin 
phases for pulsar B, which are necessary to the eclipse fitting. 
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Fig. 1 Schematic view of the double pulsar system showing the important parameters for the 
modeling of pulsar A's eclipse (dimensions and angles are not to scale). Pulsar B is located at the 
origin of the cartesian coordinate system while the projected orbital motion of pulsar A during 
its eclipse is parallel to the y axis at a constant z as seen from Earth, which is located toward 
the positive x axis. Note that since the orbital inclination is almost perfectly edge-on (|7?)). we 
can approximate the z axis to be coincident with the orbital angular momentum. The spin axis 
of pulsar B, whose spatial orientation is described by 8 and 0, is represented by the fi vector. 
The magnetic axis of pulsar B corresponds to the fi vector and makes an angle a with respect to 
fi. Finally, the absorbing region of the dipolar magnetosphere of pulsar B, truncated at radius 
-Rmag, is shown as a shaded red region. 

Fig. 2 Evolution of pulsar B's geometry as a function of time. The marginalized posterior 
probability distribution of the magnetic inclination (a), the colatitude of the spin axis (8) and 
the longitude of the spin axis (0) of pulsar B are shown from top to bottom, respectively. For 
each data point, the circle represents the median value of the posterior probability density while 
the box and the bar indicate the la and 3a confidence intervals, respectively. The gray regions 
are the 3a confidence regions derived from the joint time-dependent model fitting. Note that for 
clarity, multiple eclipses are displayed as single data points when observed over an interval of 
about a week. 

Fig. 3 Average eclipse profile of pulsar A consisting of eight eclipses observed at 820 MHz 
over a five-day period around April 11, 2007 (black line) along with a model eclipse profile 
(red dashed line). The relative pulsed flux density of pulsar A is normalized so the average 
level outside the eclipse region is unity. The resolution of each data point is ~91 ms while 1° in 
orbital phase corresponds to 24.5 s. Note that near orbital phase 0.0 the spikes are separated by 
the spin period of pulsar B. 

Fig. 4 Mass-mass diagram illustrating the present tests constraining general relativity in the 
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double pulsar system. The inset shows an expanded view of the region where the lines intersect. 
If general relativity is the correct theory of gravity, all lines should intersect at common values 
of masses. The mass ratio (R = xb/%a) an d five post-Keplerian parameters (s and r, Shapiro 
delay 'shape' and 'range'; u, periastron advance; P b , orbital period decay due to the emission 
of gravitational waves; and 7, gravitational redshift and time dilation) were reported in (|7?)) . 
Shaded regions are unphysical solutions since sini < 1, where i is the orbital inclination. In 
addition to allowing a test of the strong-field parameter ( £ -% B - J , the spin precession rate of pulsar 
B, fi B , yields a new constraint on the mass-mass diagram. 
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Parameter 


Mean 


Median 


68.2% Confidence 


99.7% Confidence 


a 


70.92° 


70.94° 


[70.49,71.31]° 


[69.68,72.13]° 


00 


130.02° 


130.02° 


[129.58, 130.44]° 


[128.79, 131.37]° 


<Pn 


51.21° 


51.20° 


[50.39, 52.03]° 


[48.80, 53.72]° 


^B 


4.77°yr- 1 


4.76°yr~ 1 


[4.12, 5.43]°yr- 1 


[2.89, 6.90]°yr- 1 



Table 1: Geometrical parameters of pulsar B derived from the eclipse model fitting. Note that 
the presented values include priors related to systematic uncertainties. The epoch of = O is 
May 2, 2006 (MJD 53857). 
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Supporting Online Material 

1 Eclipse Model 

For the eclipse model we follow the prescription of ©, section 5.2, which proposes that the 
eclipses result from synchrotron absorption of relativistic electrons trapped in the truncated 
dipolar magnetic field of pulsar B. The intensity of the transmitted radio emission from pulsar 
A is calculated from the optical depth of pulsar B's magnetosphere, r. For a given data point in 
the light curve, the optical depth is calculated as follows: 

rRmag {Bsmn\ ( x \ 
r = H [-£ ) d [-5 ' (1) 

where \i is a free parameter accounting for the characteristic optical depth of the magnetosphere 
at 820 MHz; x is the radial position along our line of sight in units of -R mag , the truncation radius 
of the dipole magnetosphere; B is the local magnetic dipolar field strength in units of -£> mag , the 
strength at the truncation radius; and k denotes the angle between our line of sight and the local 
direction of the magnetic field. Beyond the truncation radius, the optical depth is assumed to 
be zero. Pulsar B is at the origin of a coordinate system chosen so that x is along our line of 
sight, y is in the plane of the sky along the projected orbital motion and z is in the plane of the 
sky orthogonal to the x — y plane (i.e. z is parallel to the projected angular orbital momentum). 
In this coordinate system the apparent motion of pulsar A is at a fixed value z = z . The 
orientation of the pulsar B spin axis can be described by two angles, 9 the colatitude of the spin 
axis with respect to the z axis, and the longitude of the spin axis with respect to the x — z 
plane. Finally, the inclination of the magnetic axis with respect to the spin axis of pulsar B is a 
(see Fig. 1 in the main article for a schematic view of the double pulsar showing the important 
model parameters). 

Relativistic spin precession is expected to cause the spin angular momentum vector of pulsar 
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B to precess around the total angular momentum of the system, which we take to be the same 
as the orbital angular momentum since the spin contribution is negligible (~0.0001%). The 
time evolution of the pulsar spin axis in a coordinate system aligned with the orbital angular 
momentum is given by: 

S = 5 , (2) 

0so = 0so o - ^B* , (3) 

where 5 and SO are the colatitude and the longitude of the spin axis with respect to the total 
angular momentum, respectively (2). They are related to our coordinate system in the following 
way: 

cos 9 = cos(90° — i) cos 5 — sin(90° — i) sin 5 cos (f) so , (4) 

, sin5sin0 so 
sm0 = — . (5) 

ant/ 

Because we observe the system almost perfectly edge-on (i ps 90°) the z axis and the total 
angular momentum are almost perfectly coincident. Hence, we can make the approximation 
9^5, and ps so . Therefore, the evolution of the spin geometry given by equations [2] and [3] 
holds in our coordinate system as well. 

2 Eclipse Model Fitting 

Our Bayesian fitting of the model to individual eclipse data involved the six free parameters 
introduced above: \x, z , R m3ig , a, 9 and (see Fig. 1 in the main article for a schematic view of 
the model geometry). We assumed flat priors for all parameters and evaluated the joint posterior 
probability using a Markov Chain Monte Carlo (3) for the exploratory phase of the eclipse mod- 
elling. We coupled the Markov Chain Monte Carlo to simulated annealing (4) in the early stage 
of the chain to boost up convergence and explore the parameter space more efficiently. We ran 
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between three and five independent Markov Chains for each fit in order to check for consistency 
in the results. With this fitting optimization technique, we explored the full parameter sub-space 
corresponding to the geometric orientation of pulsar B's magnetosphere (i.e. a, 9 and <p) as well 
as enough of the other parameters to identify the best-fit values. 

A quantitative treatment of changes in the eclipses requires incorporation of the eclipse 
model in a framework accounting for parameter evolution. Except for <p, we do not find any sig- 
nificant secular evolution of the model parameters from their marginal posterior probability ©. 
Indeed, theoretical considerations let us presume that only should change over time due to 
relativistic spin precession. The Markov Chain Monte Carlo technique poorly samples regions 
of low probability and so we fixed \x — 2, -R mag = 1.29° (projected value in terms of orbital 
phase) and z /R mag = —0.543, their best-fit values, before making a deeper investigation of 
pulsar B's spin axis evolution. Although 9 and a are not expected to vary, we did not fix them 
as they are highly covariant with <p, as opposed to the other parameters, but also because they 
contain valuable information about pulsar B's geometry. Then, we completed our analysis of 
the eclipse evolution using a high resolution grid to map the now-reduced parameter space like- 
lihood, containing three free parameters. Following equations |2] and |3j this analysis implements 
a Bayesian joint fit that assumes a = a and 9 = 9 are constant, and that linearly varies with 
time, i.e. <p = 0o — ^b^, where fie is the precession rate. The epoch of = 0o is May 2, 2006 
(MJD 53857). Detailed results from this fit showing different combinations of one- and two- 
dimensional marginal posterior probabilities are presented in Fig. SI. Values reported in Table 
1 of the main article were derived from these one-dimensional marginal posterior probability 
functions, assuming that the integrated probabilities are unity. 
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3 Analysis of Systematics 

We investigated the importance of systematics in the eclipse modeling and concluded that two 
main effects should contribute to increasing the total uncertainty in our best-fit geometric pa- 
rameters, above the statistical value. First, we observe considerable changes in the pulse profile 
of pulsar B as shown in Fig. S2. Since the spin phases of pulsar B are input data for the mod- 
eling, losing the fiducial reference to the neutron star surface will introduce additional error in 
the fitted eclipse parameters. While we do not require a measurement of the spin phases as ac- 
curate as for timing purposes, a few percent offset translates into slightly different geometrical 
parameters. The main effect of varying spin phases is to assign earlier or later rotational phases 
that mimic a slightly faster or slower precession rate. The pulse profile evolution of pulsar B is 
likely caused by the changing viewing geometry due to relavistic spin precession. Although it 
is not clear how pulsar B's pulse profile geometry is related to its surface, we are confident that 
the technique we used to determine the spin phases yields reliable results. 

A second source of systematics arises from the choice of the eclipse region to include in the 
fit. Changes in the eclipse light curve due to relativistic spin precession are not uniform and the 
eclipse model tends to perform better toward the eclipse center than at the ingress or the egress. 
As opposed to the eclipse center, where our sight line to pulsar A goes deep inside and outside 
the magnetosphere of pulsar B as it rotates, our sight line only briefly intersects the edge of the 
magnetosphere at the beginning and the end of the eclipse. Therefore, local distortions of pulsar 
B's magnetic field or variations of the plasma density may give rise to a slight departure from our 
model. Indeed, we observe that fitting the whole eclipse does not generally provide qualitatively 
good fits. The narrow and periodic modulations in the eclipse center are very important markers 
for the geometric orientation of pulsar B but they tend to be misfitted because broader features 
in the egress region lead to larger variations of the goodness-of-fit. We find that excluding 
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the egress more accurately fits the overall light curves, without sacrificing critical information 
derived from narrow modulations, while still qualitatively reproducing the egress. Therefore, 
we chose to fit the eclipse in the range [—1.0°, 0.75°] centered around conjunction (see Fig. 3 of 
the main article for an example of a fitted eclipse light curve). 

Determining the boundaries of the region to fit is arbitrary and hence we estimated how 
much dispersion in the best- fit values is induced by other choices of limits. We compared our 
actual choice, [-1.0°, 0.75°], with the full eclipse, [-1.0°, 1.0°], the eclipse center, [-0.6°, 0.6°], 
and the extended center, [—0.7°, 0.7°], fits. In a Bayesian framework, we can easily incorporate 
the effect of systematics as priors on the model parameters. For simplicity and because the 
functional form of the systematics is poorly defined we assume Gaussian priors. Therefore, we 
can recast the posterior probability distribution of our pre-systematics analysis work, for which 
we were assuming constant priors: 

p(a,6,(f>\D)<x£(D\a,6,(j>) , (6) 

as: 

p(a,9,<j>\D) oc / Af(a-a',a a ) f N(9-9',a e ) I MU-(f>', (7 J C (Die/, 9', <j>') da' d6' d(f>' , 
Jo! Je 1 J<t>' 

(7) 
where Af(v, cr u ) is a Gaussian distribution of mean v and standard deviation o v . The likelihood, 
C (D\a, 9, 0), is defined as exp(— x^/2), with xl being the standard reduced chi-square. From 
the analysis of systematics due to the choice of the region to fit, along with the additional un- 
certainty in the spin phase of pulsar B due to the long-term pulse profile variations, we estimate 
that systematics contribute a a = 1°, ag = 1° and 0$ = 2.0°. Note that incorporating Gaussian 
priors due to systematics has the effect of convolving the three-dimensional likelihood obtained 
from the eclipse fitting, C (D\a, 9, 0), with a three-dimensional Gaussian. The analysis reported 
in this article includes these priors (this is particularly relevant for numbers quoted in Table 1 of 
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the main article and for the marginalized posterior probability distributions presented in Fig. SI 
of this Supporting Online Material). 
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Fig. SI Pulsar B's pulse profile consisting of the integrated flux at all orbital phases for 
each observation in our data set. Pulse profiles are normalized so the peak is unity and they 
are displayed with an incremental 1-unit vertical shift for clarity (the vertical axis does not 
show a linear time sequence). The vertical red dashed line marks the fiducial spin phase, which 
was determined by aligning the profiles using the first ten Fourier bins of the original 512-bin 
profile assuming that the two narrow peaks visible in the more recent data are the "edges" of 
an underlying unimodal envelope reminiscent of the profile in the earlier observations. Note 
that the observation length and radio interference contamination slightly varies from epoch to 
epoch but the overall signal-to-noise ratio decreases in the latest observations due to pulsar B 
becoming weaker. MJD 52996 corresponds to Dec. 23, 2003, and MJD 54430 to Nov. 26, 
2007. 

Fig. S2 One- and two-dimensional projections of the marginalized posterior probability 
distributions for the joint fit of the parameters' evolution. Black contours in two-dimensional 
maps are joint 1, 2, 3, 4 and 5c confidence regions, with the red color being associated with a 
higher likelihood value. The epoch of = O is May 2, 2006 (MJD 53857). Note that these 
probability distributions include priors related to systematic uncertainties (see Section 3 of this 
Supporting Online Material). 
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Movie SILj The eclipses in the double pulsar PSR J0737— 3039A/B occur when pulsar A's 
projected orbital motion, represented by a gray circle moving on a black line, passes behind its 
companion, pulsar B. Radio emission from pulsar A is absorbed via synchrotron resonance with 
the plasma trapped in the closed field lines of the truncated dipolar magnetosphere of pulsar B, 
shown as a colored dipolar structure. Since pulsar B's magnetic dipole axis is misaligned with 
respect to its spin axis (represented by a diagonal rod), the opacity along our sight line to pulsar 
A varies as a function of pulsar B's spin phase. The theoretical light curve resulting from the 
eclipse animated in the upper panel is drawn as a black curve in the bottom panel and real 
eclipse data, observed with the Green Bank Telescope in April 2007, are overlaid in red. 

Movie S2 1 Time-lapse animation displaying the evolution of pulsar B's geometry in the 
double pulsar PSR J0737— 3039A/B due to relativistic spin precession between January 2004 
and January 2029. The truncated dipolar magnetosphere of pulsar B, shown as a colored dipolar 
structure, rotates about its spin axis, pictured as a diagonal rod. The apparent orbital motion 
of pulsar A during the eclipse corresponds to the horizontal black line intersecting pulsar B's 
magnetosphere. Relativistic spin precession is similar to the wobbling of a spinning top and 
induces a motion of the spin-axis orientation around the orbital angular momentum, which is 
vertical in this movie. The theoretical light curve corresponding to the eclipse animated in the 
upper panel is drawn in the lower panel. The angle corresponds to the longitude of the spin 
axis, with 0° being the direction coincident with the line of sight. 



1 Movies, as well as the published version of this paper, can be found at 

|http: //www. sciencemag. org/cgi/content/full/sci; 321/5885/104/ DC 1| 
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